home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209b.zip / octave-2.09 / doc / octave.i07 (.txt) < prev    next >
GNU Info File  |  1997-08-20  |  50KB  |  1,104 lines

  1. This is Info file octave, produced by Makeinfo-1.64 from the input file
  2. octave.tex.
  3. START-INFO-DIR-ENTRY
  4. * Octave: (octave).     Interactive language for numerical computations.
  5. END-INFO-DIR-ENTRY
  6.    Copyright (C) 1996, 1997 John W. Eaton.
  7.    Permission is granted to make and distribute verbatim copies of this
  8. manual provided the copyright notice and this permission notice are
  9. preserved on all copies.
  10.    Permission is granted to copy and distribute modified versions of
  11. this manual under the conditions for verbatim copying, provided that
  12. the entire resulting derived work is distributed under the terms of a
  13. permission notice identical to this one.
  14.    Permission is granted to copy and distribute translations of this
  15. manual into another language, under the above conditions for modified
  16. versions.
  17. File: octave,  Node: Rearranging Matrices,  Next: Special Utility Matrices,  Prev: Finding Elements and Checking Conditions,  Up: Matrix Manipulation
  18. Rearranging Matrices
  19. ====================
  20.  - Function File:  fliplr (X)
  21.      Return a copy of X with the order of the columns reversed.  For
  22.      example,
  23.           fliplr ([1, 2; 3, 4])
  24.                =>  2  1
  25.                    4  3
  26.  - Function File:  flipud (X)
  27.      Return a copy of X with the order of the rows reversed.  For
  28.      example,
  29.           flipud ([1, 2; 3, 4])
  30.                =>  3  4
  31.                    1  2
  32.  - Function File:  rot90 (X, N)
  33.      Return a copy of X with the elements rotated counterclockwise in
  34.      90-degree increments.  The second argument is optional, and
  35.      specifies how many 90-degree rotations are to be applied (the
  36.      default value is 1).  Negative values of N rotate the matrix in a
  37.      clockwise direction.  For example,
  38.           rot90 ([1, 2; 3, 4], -1)
  39.                =>  3  1
  40.                    4  2
  41.      rotates the given matrix clockwise by 90 degrees.  The following
  42.      are all equivalent statements:
  43.           rot90 ([1, 2; 3, 4], -1)
  44.           ==
  45.           rot90 ([1, 2; 3, 4], 3)
  46.           ==
  47.           rot90 ([1, 2; 3, 4], 7)
  48.  - Function File:  reshape (A, M, N)
  49.      Return a matrix with M rows and N columns whose elements are taken
  50.      from the matrix A.  To decide how to order the elements, Octave
  51.      pretends that the elements of a matrix are stored in column-major
  52.      order (like Fortran arrays are stored).
  53.      For example,
  54.           reshape ([1, 2, 3, 4], 2, 2)
  55.                =>  1  3
  56.                    2  4
  57.      If the variable `do_fortran_indexing' is nonzero, the `reshape'
  58.      function is equivalent to
  59.           retval = zeros (m, n);
  60.           retval (:) = a;
  61.      but it is somewhat less cryptic to use `reshape' instead of the
  62.      colon operator.  Note that the total number of elements in the
  63.      original matrix must match the total number of elements in the new
  64.      matrix.
  65.  - Function File:  shift (X, B)
  66.      If X is a vector, perform a circular shift of length B of the
  67.      elements of X.
  68.      If X is a matrix, do the same for each column of X.
  69.  - Loadable Function: [S, I] = sort (X)
  70.      Return a copy of X with the elements elements arranged in
  71.      increasing order.  For matrices, `sort' orders the elements in each
  72.      column.
  73.      For example,
  74.           sort ([1, 2; 2, 3; 3, 1])
  75.                =>  1  1
  76.                    2  2
  77.                    3  3
  78.      The `sort' function may also be used to produce a matrix
  79.      containing the original row indices of the elements in the sorted
  80.      matrix.  For example,
  81.           [s, i] = sort ([1, 2; 2, 3; 3, 1])
  82.                => s = 1  1
  83.                       2  2
  84.                       3  3
  85.                => i = 1  3
  86.                       2  1
  87.                       3  2
  88.    Since the `sort' function does not allow sort keys to be specified,
  89. it can't be used to order the rows of a matrix according to the values
  90. of the elements in various columns(1) in a single call.  Using the
  91. second output, however, it is possible to sort all rows based on the
  92. values in a given column.  Here's an example that sorts the rows of a
  93. matrix based on the values in the second column.
  94.      a = [1, 2; 2, 3; 3, 1];
  95.      [s, i] = sort (a (:, 2));
  96.      a (i, :)
  97.           =>  3  1
  98.               1  2
  99.               2  3
  100.  - Function File:  tril (A, K)
  101.  - Function File:  triu (A, K)
  102.      Return a new matrix formed by extracting extract the lower (`tril')
  103.      or upper (`triu') triangular part of the matrix A, and setting all
  104.      other elements to zero.  The second argument is optional, and
  105.      specifies how many diagonals above or below the main diagonal
  106.      should also be set to zero.
  107.      The default value of K is zero, so that `triu' and `tril' normally
  108.      include the main diagonal as part of the result matrix.
  109.      If the value of K is negative, additional elements above (for
  110.      `tril') or below (for `triu') the main diagonal are also selected.
  111.      The absolute value of K must not be greater than the number of
  112.      sub- or super-diagonals.
  113.      For example,
  114.           tril (ones (3), -1)
  115.                =>  0  0  0
  116.                    1  0  0
  117.                    1  1  0
  118.      and
  119.           tril (ones (3), 1)
  120.                =>  1  1  0
  121.                    1  1  1
  122.                    1  1  1
  123.  - Function File:  vec (X)
  124.      Return the vector obtained by stacking the columns of the matrix X
  125.      one above the other.
  126.  - Function File:  vech (X)
  127.      Return the vector obtained by eliminating all supradiagonal
  128.      elements of the square matrix X and stacking the result one column
  129.      above the other.
  130.    ---------- Footnotes ----------
  131.    (1)  For example, to first sort based on the values in column 1, and
  132. then, for any values that are repeated in column 1, sort based on the
  133. values found in column 2, etc.
  134. File: octave,  Node: Special Utility Matrices,  Next: Famous Matrices,  Prev: Rearranging Matrices,  Up: Matrix Manipulation
  135. Special Utility Matrices
  136. ========================
  137.  - Built-in Function:  eye (X)
  138.  - Built-in Function:  eye (N, M)
  139.      Return an identity matrix.  If invoked with a single scalar
  140.      argument, `eye' returns a square matrix with the dimension
  141.      specified.  If you supply two scalar arguments, `eye' takes them
  142.      to be the number of rows and columns.  If given a vector with two
  143.      elements, `eye' uses the values of the elements as the number of
  144.      rows and columns, respectively.  For example,
  145.           eye (3)
  146.                =>  1  0  0
  147.                    0  1  0
  148.                    0  0  1
  149.      The following expressions all produce the same result:
  150.           eye (2)
  151.           ==
  152.           eye (2, 2)
  153.           ==
  154.           eye (size ([1, 2; 3, 4])
  155.      For compatibility with MATLAB, calling `eye' with no arguments is
  156.      equivalent to calling it with an argument of 1.
  157.  - Built-in Function:  ones (X)
  158.  - Built-in Function:  ones (N, M)
  159.      Return a matrix whose elements are all 1.  The arguments are
  160.      handled the same as the arguments for `eye'.
  161.      If you need to create a matrix whose values are all the same, you
  162.      should use an expression like
  163.           val_matrix = val * ones (n, m)
  164.  - Built-in Function:  zeros (X)
  165.  - Built-in Function:  zeros (N, M)
  166.      Return a matrix whose elements are all 0.  The arguments are
  167.      handled the same as the arguments for `eye'.
  168.  - Loadable Function:  rand (X)
  169.  - Loadable Function:  rand (N, M)
  170.  - Loadable Function:  rand (`"seed"', X)
  171.      Return a matrix with random elements uniformly distributed on the
  172.      interval (0, 1).  The arguments are handled the same as the
  173.      arguments for `eye'.  In addition, you can set the seed for the
  174.      random number generator using the form
  175.           rand ("seed", X)
  176.      where X is a scalar value.  If called as
  177.           rand ("seed")
  178.      `rand' returns the current value of the seed.
  179.  - Loadable Function:  randn (X)
  180.  - Loadable Function:  randn (N, M)
  181.  - Loadable Function:  randn (`"seed"', X)
  182.      Return a matrix with normally distributed random elements.  The
  183.      arguments are handled the same as the arguments for `eye'.  In
  184.      addition, you can set the seed for the random number generator
  185.      using the form
  186.           randn ("seed", X)
  187.      where X is a scalar value.  If called as
  188.           randn ("seed")
  189.      `randn' returns the current value of the seed.
  190.    The `rand' and `randn' functions use separate generators.  This
  191. ensures that
  192.      rand ("seed", 13);
  193.      randn ("seed", 13);
  194.      u = rand (100, 1);
  195.      n = randn (100, 1);
  196.      rand ("seed", 13);
  197.      randn ("seed", 13);
  198.      u = zeros (100, 1);
  199.      n = zeros (100, 1);
  200.      for i = 1:100
  201.        u(i) = rand ();
  202.        n(i) = randn ();
  203.      end
  204. produce equivalent results.
  205.    Normally, `rand' and `randn' obtain their initial seeds from the
  206. system clock, so that the sequence of random numbers is not the same
  207. each time you run Octave.  If you really do need for to reproduce a
  208. sequence of numbers exactly, you can set the seed to a specific value.
  209.    If it is invoked without arguments, `rand' and `randn' return a
  210. single element of a random sequence.
  211.    The `rand' and `randn' functions use Fortran code from RANLIB, a
  212. library of fortran routines for random number generation, compiled by
  213. Barry W. Brown and James Lovato of the Department of Biomathematics at
  214. The University of Texas, M.D. Anderson Cancer Center, Houston, TX 77030.
  215.  - Built-in Function:  diag (V, K)
  216.      Return a diagonal matrix with vector V on diagonal K.  The second
  217.      argument is optional.  If it is positive, the vector is placed on
  218.      the K-th super-diagonal.  If it is negative, it is placed on the
  219.      -K-th sub-diagonal.  The default value of K is 0, and the vector
  220.      is placed on the main diagonal.  For example,
  221.           diag ([1, 2, 3], 1)
  222.                =>  0  1  0  0
  223.                    0  0  2  0
  224.                    0  0  0  3
  225.                    0  0  0  0
  226.    The functions `linspace' and `logspace' make it very easy to create
  227. vectors with evenly or logarithmically spaced elements.  *Note Ranges::.
  228.  - Function File:  linspace (BASE, LIMIT, N)
  229.      Return a row vector with N linearly spaced elements between BASE
  230.      and LIMIT.  The number of elements, N, must be greater than 1.
  231.      The BASE and LIMIT are always included in the range.  If BASE is
  232.      greater than LIMIT, the elements are stored in decreasing order.
  233.      If the number of points is not specified, a value of 100 is used.
  234.      The `linspace' function always returns a row vector, regardless of
  235.      the value of `prefer_column_vectors'.
  236.  - Function File:  logspace (BASE, LIMIT, N)
  237.      Similar to `linspace' except that the values are logarithmically
  238.      spaced from 10^base to 10^limit.
  239.      If LIMIT is equal to pi, the points are between 10^base and pi,
  240.      *not* 10^base and 10^pi, in order to  be compatible with the
  241.      corresponding MATLAB function.
  242.  - Built-in Variable: treat_neg_dim_as_zero
  243.      If the value of `treat_neg_dim_as_zero' is nonzero, expressions
  244.      like
  245.           eye (-1)
  246.      produce an empty matrix (i.e., row and column dimensions are zero).
  247.      Otherwise, an error message is printed and control is returned to
  248.      the top level.  The default value is 0.
  249. File: octave,  Node: Famous Matrices,  Prev: Special Utility Matrices,  Up: Matrix Manipulation
  250. Famous Matrices
  251. ===============
  252.    The following functions return famous matrix forms.
  253.  - Function File:  hadamard (K)
  254.      Return the Hadamard matrix of order n = 2^k.
  255.  - Function File:  hankel (C, R)
  256.      Return the Hankel matrix constructed given the first column C, and
  257.      (optionally) the last row R.  If the last element of C is not the
  258.      same as the first element of R, the last element of C is used.  If
  259.      the second argument is omitted, the last row is taken to be the
  260.      same as the first column.
  261.      A Hankel matrix formed from an m-vector C, and an n-vector R, has
  262.      the elements
  263.           H (i, j) = c (i+j-1),  i+j-1 <= m;
  264.           H (i, j) = r (i+j-m),  otherwise
  265.  - Function File:  hilb (N)
  266.      Return the Hilbert matrix of order N.  The i, j element of a
  267.      Hilbert matrix is defined as
  268.           H (i, j) = 1 / (i + j - 1)
  269.  - Function File:  invhilb (N)
  270.      Return the inverse of a Hilbert matrix of order N.  This is exact.
  271.      Compare with the numerical calculation of `inverse (hilb (n))',
  272.      which suffers from the ill-conditioning of the Hilbert matrix, and
  273.      the finite precision of your computer's floating point arithmetic.
  274.  - Function File:  toeplitz (C, R)
  275.      Return the Toeplitz matrix constructed given the first column C,
  276.      and (optionally) the first row R.  If the first element of C is
  277.      not the same as the first element of R, the first element of C is
  278.      used.  If the second argument is omitted, the first row is taken
  279.      to be the same as the first column.
  280.      A square Toeplitz matrix has the form
  281.           c(0)  r(1)   r(2)  ...  r(n)
  282.           c(1)  c(0)   r(1)      r(n-1)
  283.           c(2)  c(1)   c(0)      r(n-2)
  284.            .                       .
  285.            .                       .
  286.            .                       .
  287.           
  288.           c(n) c(n-1) c(n-2) ...  c(0)
  289.  - Function File:  vander (C)
  290.      Return the Vandermonde matrix whose next to last column is C.
  291.      A Vandermonde matrix has the form
  292.           c(0)^n ... c(0)^2  c(0)  1
  293.           c(1)^n ... c(1)^2  c(1)  1
  294.            .           .      .    .
  295.            .           .      .    .
  296.            .           .      .    .
  297.           
  298.           c(n)^n ... c(n)^2  c(n)  1
  299. File: octave,  Node: Arithmetic,  Next: Linear Algebra,  Prev: Matrix Manipulation,  Up: Top
  300. Arithmetic
  301. **********
  302.    Unless otherwise noted, all of the functions described in this
  303. chapter will work for real and complex scalar or matrix arguments.
  304. * Menu:
  305. * Utility Functions::
  306. * Complex Arithmetic::
  307. * Trigonometry::
  308. * Sums and Products::
  309. * Special Functions::
  310. * Mathematical Constants::
  311. File: octave,  Node: Utility Functions,  Next: Complex Arithmetic,  Prev: Arithmetic,  Up: Arithmetic
  312. Utility Functions
  313. =================
  314.    The following functions are available for working with complex
  315. numbers.  Each expects a single argument.  They are called "mapping
  316. functions" because when given a matrix argument, they apply the given
  317. function to each element of the matrix.
  318.  - Mapping Function:  ceil (X)
  319.      Return the smallest integer not less than X.  If X is complex,
  320.      return `ceil (real (X)) + ceil (imag (X)) * I'.
  321.  - Mapping Function:  exp (X)
  322.      Compute the exponential of X.  To compute the matrix exponential,
  323.      see *Note Linear Algebra::.
  324.  - Mapping Function:  fix (X)
  325.      Truncate X toward zero.  If X is complex, return `fix (real (X)) +
  326.      fix (imag (X)) * I'.
  327.  - Mapping Function:  floor (X)
  328.      Return the largest integer not greater than X.  If X is complex,
  329.      return `floor (real (X)) + floor (imag (X)) * I'.
  330.  - Mapping Function:  gcd (X, `...')
  331.      Compute the greatest common divisor of the elements of X, or the
  332.      list of all the arguments.  For example,
  333.           gcd (a1, ..., ak)
  334.      is the same as
  335.           gcd ([a1, ..., ak])
  336.      An optional second return value, V contains an integer vector such
  337.      that
  338.           g = v(1) * a(k) + ... + v(k) * a(k)
  339.  - Mapping Function:  lcm (X, `...')
  340.      Compute the least common multiple of the elements elements of X, or
  341.      the list of all the arguments.  For example,
  342.           lcm (a1, ..., ak)
  343.      is the same as
  344.           lcm ([a1, ..., ak]).
  345.  - Mapping Function:  log (X)
  346.      Compute the natural logarithm of X.  To compute the matrix
  347.      logarithm, see *Note Linear Algebra::.
  348.  - Mapping Function:  log10 (X)
  349.      Compute the base-10 logarithm of X.
  350.  - Mapping Function: Y = log2 (X)
  351.  - Mapping Function: [F, E] log2 (X)
  352.      Compute the base-2 logarithm of X.  With two outputs, returns F
  353.      and E such that  1/2 <= abs(f) < 1 and x = f * 2^e.
  354.  - Loadable Function:  max (X)
  355.      For a vector argument, return the maximum value.  For a matrix
  356.      argument, return the maximum value from each column, as a row
  357.      vector.  Thus,
  358.           max (max (X))
  359.      returns the largest element of X.
  360.      For complex arguments, the magnitude of the elements are used for
  361.      comparison.
  362.  - Loadable Function:  min (X)
  363.      Like `max', but return the minimum value.
  364.  - Function File:  nextpow2 (X)
  365.      If X is a scalar, returns the first integer N such that  2^n >=
  366.      abs (x).
  367.      If X is a vector, return `nextpow2 (length (X))'.
  368.  - Mapping Function:  pow2 (X)
  369.  - Mapping Function:  pow2 (F, E)
  370.      With one argument, computes  2 .^ x for each element of X.  With
  371.      two arguments, returns  f .* (2 .^ e).
  372.  - Mapping Function:  rem (X, Y)
  373.      Return the remainder of `X / Y', computed using the expression
  374.           x - y .* fix (x ./ y)
  375.      An error message is printed if the dimensions of the arguments do
  376.      not agree, or if either of the arguments is complex.
  377.  - Mapping Function:  round (X)
  378.      Return the integer nearest to X.  If X is complex, return `round
  379.      (real (X)) + round (imag (X)) * I'.
  380.  - Mapping Function:  sign (X)
  381.      Compute the "signum" function, which is defined as
  382.                      -1, x < 0;
  383.           sign (x) =  0, x = 0;
  384.                       1, x > 0.
  385.      For complex arguments, `sign' returns `x ./ abs (X)'.
  386.  - Mapping Function:  sqrt (X)
  387.      Compute the square root of X.  If X is negative, a complex result
  388.      is returned.  To compute the matrix square root, see *Note Linear
  389.      Algebra::.
  390.  - Mapping Function:  xor (X, Y)
  391.      Return the `exclusive or' of the entries of X and Y.  For boolean
  392.      expressions X and Y, `xor (X, Y)' is true if and only if X or Y is
  393.      true, but not if both X and Y are true.
  394. File: octave,  Node: Complex Arithmetic,  Next: Trigonometry,  Prev: Utility Functions,  Up: Arithmetic
  395. Complex Arithmetic
  396. ==================
  397.    The following functions are available for working with complex
  398. numbers.  Each expects a single argument.  Given a matrix they work on
  399. an element by element basis.  In the descriptions of the following
  400. functions, Z is the complex number X + IY, where I is defined as `sqrt
  401. (-1)'.
  402.  - Mapping Function:  abs (Z)
  403.      Compute the magnitude of Z, defined as |Z| = `sqrt (x^2 + y^2)'.
  404.      For example,
  405.           abs (3 + 4i)
  406.                => 5
  407.  - Mapping Function:  arg (Z)
  408.  - Mapping Function:  angle (Z)
  409.      Compute the argument of Z, defined as THETA = `atan (Y/X)'.
  410.      in radians.
  411.      For example,
  412.           arg (3 + 4i)
  413.                => 0.92730
  414.  - Mapping Function:  conj (Z)
  415.      Return the complex conjugate of Z, defined as `conj (Z)' = X - IY.
  416.  - Mapping Function:  imag (Z)
  417.      Return the imaginary part of Z as a real number.
  418.  - Mapping Function:  real (Z)
  419.      Return the real part of Z.
  420. File: octave,  Node: Trigonometry,  Next: Sums and Products,  Prev: Complex Arithmetic,  Up: Arithmetic
  421. Trigonometry
  422. ============
  423.    Octave provides the following trigonometric functions:
  424.  - Mapping Function:  sin (Z)
  425.  - Mapping Function:  cos (Z)
  426.  - Mapping Function:  tan (Z)
  427.  - Mapping Function:  sec (Z)
  428.  - Mapping Function:  csc (Z)
  429.  - Mapping Function:  cot (Z)
  430.      The ordinary trigonometric functions.
  431.  - Mapping Function:  asin (Z)
  432.  - Mapping Function:  acos (Z)
  433.  - Mapping Function:  atan (Z)
  434.  - Mapping Function:  asec (Z)
  435.  - Mapping Function:  acsc (Z)
  436.  - Mapping Function:  acot (Z)
  437.      The ordinary inverse trigonometric functions.
  438.  - Mapping Function:  sinh (Z)
  439.  - Mapping Function:  cosh (Z)
  440.  - Mapping Function:  tanh (Z)
  441.  - Mapping Function:  sech (Z)
  442.  - Mapping Function:  csch (Z)
  443.  - Mapping Function:  coth (Z)
  444.      Hyperbolic trigonometric functions.
  445.  - Mapping Function:  asinh (Z)
  446.  - Mapping Function:  acosh (Z)
  447.  - Mapping Function:  atanh (Z)
  448.  - Mapping Function:  asech (Z)
  449.  - Mapping Function:  acsch (Z)
  450.  - Mapping Function:  acoth (Z)
  451.      Inverse hyperbolic trigonometric functions.
  452.    Each of these functions expect a single argument.  For matrix
  453. arguments, they work on an element by element basis.  For example,
  454.      sin ([1, 2; 3, 4])
  455.           =>  0.84147   0.90930
  456.               0.14112  -0.75680
  457.  - Mapping Function:  atan2 (Y, X)
  458.      Return the arctangent of Y/X.  The signs of the arguments are used
  459.      to determine the quadrant of the result, which is in the range
  460.      `pi' to -`pi'.
  461. File: octave,  Node: Sums and Products,  Next: Special Functions,  Prev: Trigonometry,  Up: Arithmetic
  462. Sums and Products
  463. =================
  464.  - Built-in Function:  sum (X)
  465.      For a vector argument, return the sum of all the elements.  For a
  466.      matrix argument, return the sum of the elements in each column, as
  467.      a row vector.  The sum of an empty matrix is 0 if it has no
  468.      columns, or a vector of zeros if it has no rows (*note Empty
  469.      Matrices::.).
  470.  - Built-in Function:  prod (X)
  471.      For a vector argument, return the product of all the elements.
  472.      For a matrix argument, return the product of the elements in each
  473.      column, as a row vector.  The product of an empty matrix is 1 if
  474.      it has no columns, or a vector of ones if it has no rows (*note
  475.      Empty Matrices::.).
  476.  - Built-in Function:  cumsum (X)
  477.      Return the cumulative sum of each column of X.  For example,
  478.           cumsum ([1, 2; 3, 4])
  479.                =>  1  2
  480.                    4  6
  481.  - Built-in Function:  cumprod (X)
  482.      Return the cumulative product of each column of X.  For example,
  483.           cumprod ([1, 2; 3, 4])
  484.                =>  1  2
  485.                    3  8
  486.  - Built-in Function:  sumsq (X)
  487.      For a vector argument, return the sum of the squares of all the
  488.      elements.  For a matrix argument, return the sum of the squares of
  489.      the elements in each column, as a row vector.
  490. File: octave,  Node: Special Functions,  Next: Mathematical Constants,  Prev: Sums and Products,  Up: Arithmetic
  491. Special Functions
  492. =================
  493.  - Mapping Function:  beta (A, B)
  494.      Return the Beta function,
  495.           beta (a, b) = gamma (a) * gamma (b) / gamma (a + b).
  496.  - Mapping Function:  betai (A, B, X)
  497.      Return the incomplete Beta function,
  498.                                               x
  499.                                              /
  500.           betai (a, b, x) = beta (a, b)^(-1) | t^(a-1) (1-t)^(b-1) dt.
  501.                                              /
  502.                                           t=0
  503.      If x has more than one component, both A and B must be scalars.
  504.      If X is a scalar, A and B must be of compatible dimensions.
  505.  - Mapping Function:  bincoeff (N, K)
  506.      Return the binomial coefficient of N and K, defined as
  507.            /   \
  508.            | n |    n (n-1) (n-2) ... (n-k+1)
  509.            |   |  = -------------------------
  510.            | k |               k!
  511.            \   /
  512.      For example,
  513.           bincoeff (5, 2)
  514.                => 10
  515.  - Mapping Function:  erf (Z)
  516.      Computes the error function,
  517.                                    z
  518.                                   /
  519.           erf (z) = (2/sqrt (pi)) | e^(-t^2) dt
  520.                                   /
  521.                                t=0
  522.  - Mapping Function:  erfc (Z)
  523.      Computes the complementary error function, `1 - erf (Z)'.
  524.  - Mapping Function:  erfinv (Z)
  525.      Computes the inverse of the error function,
  526.  - Mapping Function:  gamma (Z)
  527.      Computes the Gamma function,
  528.                       infinity
  529.                       /
  530.           gamma (z) = | t^(z-1) exp (-t) dt.
  531.                       /
  532.                    t=0
  533.  - Mapping Function:  gammai (A, X)
  534.      Computes the incomplete gamma function,
  535.                                         x
  536.                               1        /
  537.           gammai (a, x) = ---------    | exp (-t) t^(a-1) dt
  538.                           gamma (a)    /
  539.                                     t=0
  540.      If A is scalar, then `gammai (A, X)' is returned for each element
  541.      of X and vice versa.
  542.      If neither A nor X is scalar, the sizes of A and X must agree, and
  543.      GAMMAI is applied element-by-element.
  544.  - Mapping Function:  lgamma (A, X)
  545.  - Mapping Function:  gammaln (A, X)
  546.      Return the natural logarithm of the gamma function.
  547.  - Function File:  cross (X, Y)
  548.      Computes the vector cross product of the two 3-dimensional vectors
  549.      X and Y.  For example,
  550.           cross ([1,1,0], [0,1,1])
  551.                => [ 1; -1; 1 ]
  552.  - Function File:  commutation_matrix (M, N)
  553.      Return the commutation matrix  K(m,n)  which is the unique  M*N by
  554.      M*N  matrix such that  K(M,N) * vec (A) = vec (A')  for all  M by N
  555.      matrices  A.
  556.      If only one argument M is given,  K(m,m)  is returned.
  557.      See Magnus and Neudecker (1988), Matrix differential calculus with
  558.      applications in statistics and econometrics.
  559.  - Function File:  duplication_matrix (N)
  560.      Return the duplication matrix  D_N  which is the unique  N^2 by
  561.      N*(N+1)/2  matrix such that  D_N \cdot vech (A) = vec (A)  for all
  562.      symmetric  N by N  matrices  A.
  563.      See Magnus and Neudecker (1988), Matrix differential calculus with
  564.      applications in statistics and econometrics.
  565. File: octave,  Node: Mathematical Constants,  Prev: Special Functions,  Up: Arithmetic
  566. Mathematical Constants
  567. ======================
  568.  - Built-in Variable: I
  569.  - Built-in Variable: J
  570.  - Built-in Variable: i
  571.  - Built-in Variable: j
  572.      A pure imaginary number, defined as   `sqrt (-1)'.  The `I' and
  573.      `J' forms are true constants, and cannot be modified.  The `i' and
  574.      `j' forms are like ordinary variables, and may be used for other
  575.      purposes.  However, unlike other variables, they once again assume
  576.      their special predefined values if they are cleared *Note Status
  577.      of Variables::.
  578.  - Built-in Variable: Inf
  579.  - Built-in Variable: inf
  580.      Infinity.  This is the result of an operation like 1/0, or an
  581.      operation that results in a floating point overflow.
  582.  - Built-in Variable: NaN
  583.  - Built-in Variable: nan
  584.      Not a number.  This is the result of an operation like 0/0, or
  585.      `Inf - Inf', or any operation with a NaN.
  586.  - Built-in Variable: pi
  587.      The ratio of the circumference of a circle to its diameter.
  588.      Internally, `pi' is computed as `4.0 * atan (1.0)'.
  589.  - Built-in Variable: e
  590.      The base of natural logarithms.  The constant  E  satisfies the
  591.      equation  `log' (E) = 1.
  592.  - Built-in Variable: eps
  593.      The machine precision.  More precisely, `eps' is the largest
  594.      relative spacing between any two adjacent numbers in the machine's
  595.      floating point system.  This number is obviously system-dependent.
  596.      On machines that support 64 bit IEEE floating point arithmetic,
  597.      `eps' is approximately  2.2204e-16.
  598.  - Built-in Variable: realmax
  599.      The largest floating point number that is representable.  The
  600.      actual value is system-dependent.  On machines that support 64 bit
  601.      IEEE floating point arithmetic, `realmax' is approximately
  602.      1.7977e+308
  603.  - Built-in Variable: realmin
  604.      The smallest floating point number that is representable.  The
  605.      actual value is system-dependent.  On machines that support 64 bit
  606.      IEEE floating point arithmetic, `realmin' is approximately
  607.      2.2251e-308
  608. File: octave,  Node: Linear Algebra,  Next: Nonlinear Equations,  Prev: Arithmetic,  Up: Top
  609. Linear Algebra
  610. **************
  611.    This chapter documents the linear algebra functions of Octave.
  612. Reference material for many of these functions may be found in Golub
  613. and Van Loan, `Matrix Computations, 2nd Ed.', Johns Hopkins, 1989, and
  614. in `LAPACK Users' Guide', SIAM, 1992.
  615. * Menu:
  616. * Basic Matrix Functions::
  617. * Matrix Factorizations::
  618. * Functions of a Matrix::
  619. File: octave,  Node: Basic Matrix Functions,  Next: Matrix Factorizations,  Prev: Linear Algebra,  Up: Linear Algebra
  620. Basic Matrix Functions
  621. ======================
  622.  - Loadable Function: AA = balance (A, OPT)
  623.  - Loadable Function: [DD, AA] = balance (A, OPT)
  624.  - Loadable Function: [CC, DD, AA, BB] = balance (A, B, OPT)
  625.      `[dd, aa] = balance (a)' returns `aa = dd \ a * dd'.  `aa' is a
  626.      matrix whose row and column norms are roughly equal in magnitude,
  627.      and `dd' = `p * d', where `p' is a permutation matrix and `d' is a
  628.      diagonal matrix of powers of two.  This allows the equilibration
  629.      to be computed without roundoff.  Results of eigenvalue
  630.      calculation are typically improved by balancing first.
  631.      `[cc, dd, aa, bb] = balance (a, b)' returns `aa = cc*a*dd' and `bb
  632.      = cc*b*dd)', where `aa' and `bb' have non-zero elements of
  633.      approximately the same magnitude and `cc' and `dd' are permuted
  634.      diagonal matrices as in `dd' for the algebraic eigenvalue problem.
  635.      The eigenvalue balancing option `opt' is selected as follows:
  636.     `"N"', `"n"'
  637.           No balancing; arguments copied, transformation(s) set to
  638.           identity.
  639.     `"P"', `"p"'
  640.           Permute argument(s) to isolate eigenvalues where possible.
  641.     `"S"', `"s"'
  642.           Scale to improve accuracy of computed eigenvalues.
  643.     `"B"', `"b"'
  644.           Permute and scale, in that order. Rows/columns of a (and b)
  645.           that are isolated by permutation are not scaled.  This is the
  646.           default behavior.
  647.      Algebraic eigenvalue balancing uses standard LAPACK routines.
  648.      Generalized eigenvalue problem balancing uses Ward's algorithm
  649.      (SIAM Journal on Scientific and Statistical Computing, 1981).
  650.  - :  cond (A)
  651.      Compute the (two-norm) condition number of a matrix. `cond (a)' is
  652.      defined as `norm (a) * norm (inv (a))', and is computed via a
  653.      singular value decomposition.
  654.  - Loadable Function:  det (A)
  655.      Compute the determinant of A using LINPACK.
  656.  - Loadable Function: LAMBDA = eig (A)
  657.  - Loadable Function: [V, LAMBDA] = eig (A)
  658.      The eigenvalues (and eigenvectors) of a matrix are computed in a
  659.      several step process which begins with a Hessenberg decomposition,
  660.      followed by a Schur decomposition, from which the eigenvalues are
  661.      apparent.  The eigenvectors, when desired, are computed by further
  662.      manipulations of the Schur decomposition.
  663.  - Loadable Function: G = givens (X, Y)
  664.  - Loadable Function: [C, S] = givens (X, Y)
  665.      Return a 2 by 2 orthogonal matrix `G = [C S; -S' C]' such that `G
  666.      [X; Y] = [*; 0]' with X and Y scalars.
  667.      For example,
  668.           givens (1, 1)
  669.                =>   0.70711   0.70711
  670.                    -0.70711   0.70711
  671.  - Loadable Function:  inv (A)
  672.  - Loadable Function:  inverse (A)
  673.      Compute the inverse of the square matrix A.
  674.  - Function File:  norm (A, P)
  675.      Compute the p-norm of the matrix A.  If the second argument is
  676.      missing, `p = 2' is assumed.
  677.      If A is a matrix:
  678.     P = `1'
  679.           1-norm, the largest column sum of A.
  680.     P = `2'
  681.           Largest singular value of A.
  682.     P = `Inf'
  683.           Infinity norm, the largest row sum of A.
  684.     P = `"fro"'
  685.           Frobenius norm of A, `sqrt (sum (diag (A' * A)))'.
  686.      If A is a vector or a scalar:
  687.     P = `Inf'
  688.           `max (abs (A))'.
  689.     P = `-Inf'
  690.           `min (abs (A))'.
  691.     other
  692.           p-norm of A, `(sum (abs (A) .^ P)) ^ (1/P)'.
  693.  - Function File:  null (A, TOL)
  694.      Return an orthonormal basis of the null space of A.
  695.      The dimension of the null space is taken as the number of singular
  696.      values of A not greater than TOL.  If the argument TOL is missing,
  697.      it is computed as
  698.           max (size (A)) * max (svd (A)) * eps
  699.  - Function File:  orth (A, TOL)
  700.      Return an orthonormal basis of the range space of A.
  701.      The dimension of the range space is taken as the number of singular
  702.      values of A greater than TOL.  If the argument TOL is missing, it
  703.      is computed as
  704.           max (size (A)) * max (svd (A)) * eps
  705.  - Function File:  pinv (X, TOL)
  706.      Return the pseudoinverse of X.  Singular values less than TOL are
  707.      ignored.
  708.      If the second argument is omitted, it is assumed that
  709.           tol = max (size (X)) * sigma_max (X) * eps,
  710.      where `sigma_max (X)' is the maximal singular value of X.
  711.  - Function File:  rank (A, TOL)
  712.      Compute the rank of A, using the singular value decomposition.
  713.      The rank is taken to be the number  of singular values of A that
  714.      are greater than the specified tolerance TOL.  If the second
  715.      argument is omitted, it is taken to be
  716.           tol = max (size (A)) * sigma (1) * eps;
  717.      where `eps' is machine precision and `sigma' is the largest
  718.      singular value of A.
  719.  - Function File:  trace (A)
  720.      Compute the trace of A, `sum (diag (A))'.
  721. File: octave,  Node: Matrix Factorizations,  Next: Functions of a Matrix,  Prev: Basic Matrix Functions,  Up: Linear Algebra
  722. Matrix Factorizations
  723. =====================
  724.  - Loadable Function:  chol (A)
  725.      Compute the Cholesky factor, R, of the symmetric positive definite
  726.      matrix A, where
  727.           r' * r = a.
  728.  - Loadable Function: H = hess (A)
  729.  - Loadable Function: [P, H] = hess (A)
  730.      Compute the Hessenberg decomposition of the matrix A.
  731.      The Hessenberg decomposition is usually used as the first step in
  732.      an eigenvalue computation, but has other applications as well (see
  733.      Golub, Nash, and Van Loan, IEEE Transactions on Automatic Control,
  734.      1979.  The Hessenberg decomposition is `p * h * p' = a' where `p'
  735.      is a square unitary matrix (`p' * p = I', using complex-conjugate
  736.      transposition) and `h' is upper Hessenberg (`i >= j+1 => h (i, j)
  737.      = 0').
  738.  - Loadable Function: [L, U, P] = lu (A)
  739.      Compute the LU decomposition of A, using subroutines from LAPACK.
  740.      The result is returned in a permuted form, according to the
  741.      optional return value P.  For example, given the matrix `a = [1,
  742.      2; 3, 4]',
  743.           [l, u, p] = lu (a)
  744.      returns
  745.           l =
  746.           
  747.             1.00000  0.00000
  748.             0.33333  1.00000
  749.           
  750.           u =
  751.           
  752.             3.00000  4.00000
  753.             0.00000  0.66667
  754.           
  755.           p =
  756.           
  757.             0  1
  758.             1  0
  759.  - Loadable Function: [Q, R, P] = qr (A)
  760.      Compute the QR factorization of A, using standard LAPACK
  761.      subroutines.  For example, given the matrix `a = [1, 2; 3, 4]',
  762.           [q, r] = qr (a)
  763.      returns
  764.           q =
  765.           
  766.             -0.31623  -0.94868
  767.             -0.94868   0.31623
  768.           
  769.           r =
  770.           
  771.             -3.16228  -4.42719
  772.              0.00000  -0.63246
  773.      The `qr' factorization has applications in the solution of least
  774.      squares problems
  775.           `min norm(A x - b)'
  776.      for overdetermined systems of equations (i.e., `a'  is a tall,
  777.      thin matrix).  The QR factorization is `q * r = a' where `q' is an
  778.      orthogonal matrix and `r' is upper triangular.
  779.      The permuted QR factorization `[Q, R, P] = qr (A)' forms the QR
  780.      factorization such that the diagonal entries of `r' are decreasing
  781.      in magnitude order.  For example, given the matrix `a = [1, 2; 3,
  782.      4]',
  783.           [q, r, pi] = qr(a)
  784.      returns
  785.           q =
  786.           
  787.             -0.44721  -0.89443
  788.             -0.89443   0.44721
  789.           
  790.           r =
  791.           
  792.             -4.47214  -3.13050
  793.              0.00000   0.44721
  794.           
  795.           p =
  796.           
  797.              0  1
  798.              1  0
  799.      The permuted `qr' factorization `[q, r, p] = qr (a)' factorization
  800.      allows the construction of an orthogonal basis of `span (a)'.
  801.  - Loadable Function: S = schur (A)
  802.  - Loadable Function: [U, S] = schur (A, OPT)
  803.      The Schur decomposition is used to compute eigenvalues of a square
  804.      matrix, and has applications in the solution of algebraic Riccati
  805.      equations in control (see `are' and `dare').  `schur' always
  806.      returns `s = u' * a * u' where `u'  is a unitary matrix (`u'* u'
  807.      is identity) and `s' is upper triangular.  The eigenvalues of `a'
  808.      (and `s') are the diagonal elements of `s' If the matrix `a' is
  809.      real, then the real Schur decomposition is computed, in which the
  810.      matrix `u' is orthogonal and `s' is block upper triangular with
  811.      blocks of size at most `2 x 2' blocks along the diagonal.  The
  812.      diagonal elements of `s' (or the eigenvalues of the `2 x 2'
  813.      blocks, when appropriate) are the eigenvalues of `a' and `s'.
  814.      The eigenvalues are optionally ordered along the diagonal
  815.      according to the value of `opt'.  `opt = "a"' indicates that all
  816.      eigenvalues with negative real parts should be moved to the leading
  817.      block of `s' (used in `are'), `opt = "d"' indicates that all
  818.      eigenvalues with magnitude less than one should be moved to the
  819.      leading block of `s' (used in `dare'), and `opt = "u"', the
  820.      default, indicates that no ordering of eigenvalues should occur.
  821.      The leading `k' columns of `u' always span the `a'-invariant
  822.      subspace corresponding to the `k' leading eigenvalues of `s'.
  823.  - Loadable Function: S = svd (A)
  824.  - Loadable Function: [U, S, V] = svd (A)
  825.      Compute the singular value decomposition of A
  826.           a = u * sigma * v'
  827.      The function `svd' normally returns the vector of singular values.
  828.      If asked for three return values, it computes U, S, and V.  For
  829.      example,
  830.           svd (hilb (3))
  831.      returns
  832.           ans =
  833.           
  834.             1.4083189
  835.             0.1223271
  836.             0.0026873
  837.      and
  838.           [u, s, v] = svd (hilb (3))
  839.      returns
  840.           u =
  841.           
  842.             -0.82704   0.54745   0.12766
  843.             -0.45986  -0.52829  -0.71375
  844.             -0.32330  -0.64901   0.68867
  845.           
  846.           s =
  847.           
  848.             1.40832  0.00000  0.00000
  849.             0.00000  0.12233  0.00000
  850.             0.00000  0.00000  0.00269
  851.           
  852.           v =
  853.           
  854.             -0.82704   0.54745   0.12766
  855.             -0.45986  -0.52829  -0.71375
  856.             -0.32330  -0.64901   0.68867
  857.      If given a second argument, `svd' returns an economy-sized
  858.      decomposition, eliminating the unnecessary rows or columns of U or
  859.      V.
  860. File: octave,  Node: Functions of a Matrix,  Prev: Matrix Factorizations,  Up: Linear Algebra
  861. Functions of a Matrix
  862. =====================
  863.  - Loadable Function:  expm (A)
  864.      Return the exponential of a matrix, defined as the infinite Taylor
  865.      series
  866.           expm(a) = I + a + a^2/2! + a^3/3! + ...
  867.      The Taylor series is *not* the way to compute the matrix
  868.      exponential; see Moler and Van Loan, `Nineteen Dubious Ways to
  869.      Compute the Exponential of a Matrix', SIAM Review, 1978.  This
  870.      routine uses Ward's diagonal Pade' approximation method with three
  871.      step preconditioning (SIAM Journal on Numerical Analysis, 1977).
  872.      Diagonal Pade'  approximations are rational polynomials of matrices
  873.                -1
  874.           D (a)   N (a)
  875.      whose Taylor series matches the first `2q+1' terms of the Taylor
  876.      series above; direct evaluation of the Taylor series (with the
  877.      same preconditioning steps) may be desirable in lieu of the Pade'
  878.      approximation when `Dq(a)' is ill-conditioned.
  879.  - Loadable Function:  logm (A)
  880.      Compute the matrix logarithm of the square matrix A.  Note that
  881.      this is currently implemented in terms of an eigenvalue expansion
  882.      and needs to be improved to be more robust.
  883.  - Loadable Function:  sqrtm (A)
  884.      Compute the matrix square root of the square matrix A.  Note that
  885.      this is currently implemented in terms of an eigenvalue expansion
  886.      and needs to be improved to be more robust.
  887.  - Function File:  kron (A, B)
  888.      Form the kronecker product of two matrices, defined block by block
  889.      as
  890.           x = [a(i, j) b]
  891.      For example,
  892.           kron (1:4, ones (3, 1))
  893.                =>  1  2  3  4
  894.                    1  2  3  4
  895.                    1  2  3  4
  896.  - Function File: [AA, BB, Q, Z] = qzhess (A, B)
  897.      Compute the Hessenberg-triangular decomposition of the matrix
  898.      pencil `(A, B)', returning `AA = Q * A * Z', `BB = Q * B * Z',
  899.      with Q and Z orthogonal.  For example,
  900.           [aa, bb, q, z] = qzhess ([1, 2; 3, 4], [5, 6; 7, 8])
  901.                => aa = [ -3.02244, -4.41741;  0.92998,  0.69749 ]
  902.                => bb = [ -8.60233, -9.99730;  0.00000, -0.23250 ]
  903.                =>  q = [ -0.58124, -0.81373; -0.81373,  0.58124 ]
  904.                =>  z = [ 1, 0; 0, 1 ]
  905.      The Hessenberg-triangular decomposition is the first step in Moler
  906.      and Stewart's QZ decomposition algorithm.
  907.      Algorithm taken from Golub and Van Loan, `Matrix Computations, 2nd
  908.      edition'.
  909.  - Loadable Function:  qzval (A, B)
  910.      Compute generalized eigenvalues of the matrix pencil `A - lambda
  911.      B'.
  912.      The arguments A and B must be real matrices.
  913.  - Loadable Function: X = syl (A, B, C)
  914.      Solve the Sylvester equation
  915.           A X + X B + C = 0
  916.      using standard LAPACK subroutines.  For example,
  917.           syl ([1, 2; 3, 4], [5, 6; 7, 8], [9, 10; 11, 12])
  918.                => [ -0.50000, -0.66667; -0.66667, -0.50000 ]
  919. File: octave,  Node: Nonlinear Equations,  Next: Quadrature,  Prev: Linear Algebra,  Up: Top
  920. Nonlinear Equations
  921. *******************
  922.    Octave can solve sets of nonlinear equations of the form
  923.      F (x) = 0
  924. using the function `fsolve', which is based on the MINPACK subroutine
  925. `hybrd'.
  926.  - Loadable Function: [X, INFO] = fsolve (FCN, X0)
  927.      Given FCN, the name of a function of the form `f (X)' and an
  928.      initial starting point X0, `fsolve' solves the set of equations
  929.      such that `f(X) == 0'.
  930.  - Loadable Function:  fsolve_options (OPT, VAL)
  931.      When called with two arguments, this function allows you set
  932.      options parameters for the function `fsolve'.  Given one argument,
  933.      `fsolve_options' returns the value of the corresponding option.  If
  934.      no arguments are supplied, the names of all the available options
  935.      and their current values are displayed.
  936.    Here is a complete example.  To solve the set of equations
  937.      -2x^2 + 3xy   + 4 sin(y) = 6
  938.       3x^2 - 2xy^2 + 3 cos(x) = -4
  939. you first need to write a function to compute the value of the given
  940. function.  For example:
  941.      function y = f (x)
  942.        y(1) = -2*x(1)^2 + 3*x(1)*x(2)   + 4*sin(x(2)) - 6;
  943.        y(2) =  3*x(1)^2 - 2*x(1)*x(2)^2 + 3*cos(x(1)) + 4;
  944.      endfunction
  945.    Then, call `fsolve' with a specified initial condition to find the
  946. roots of the system of equations.  For example, given the function `f'
  947. defined above,
  948.      [x, info] = fsolve ("f", [1; 2])
  949. results in the solution
  950.      x =
  951.      
  952.        0.57983
  953.        2.54621
  954.      
  955.      info = 1
  956.    A value of `info = 1' indicates that the solution has converged.
  957.    The function `perror' may be used to print English messages
  958. corresponding to the numeric error codes.  For example,
  959.      perror ("fsolve", 1)
  960.           -| solution converged to requested tolerance
  961. File: octave,  Node: Quadrature,  Next: Differential Equations,  Prev: Nonlinear Equations,  Up: Top
  962. Quadrature
  963. **********
  964. * Menu:
  965. * Functions of One Variable::
  966. * Orthogonal Collocation::
  967. File: octave,  Node: Functions of One Variable,  Next: Orthogonal Collocation,  Prev: Quadrature,  Up: Quadrature
  968. Functions of One Variable
  969. =========================
  970.  - Loadable Function: [V, IER, NFUN, ERR] = quad (F, A, B, TOL, SING)
  971.      Integrate a nonlinear function of one variable using Quadpack.
  972.      The first argument is the name of the  function to call to compute
  973.      the value of the integrand.  It must have the form
  974.           y = f (x)
  975.      where Y and X are scalars.
  976.      The second and third arguments are limits of integration.  Either
  977.      or both may be infinite.
  978.      The optional argument TOL is a vector that specifies the desired
  979.      accuracy of the result.  The first element of the vector is the
  980.      desired absolute tolerance, and the second element is the desired
  981.      relative tolerance.  To choose a relative test only, set the
  982.      absolute tolerance to zero.  To choose an absolute test only, set
  983.      the relative tolerance to zero.
  984.      The optional argument SING is a vector of values at which the
  985.      integrand is known to be singular.
  986.      The result of the integration is returned in V and IER contains an
  987.      integer error code (0 indicates a successful integration).  The
  988.      value of NFUN indicates how many function evaluations were
  989.      required, and ERR contains an estimate of the error in the
  990.      solution.
  991.  - Loadable Function:  quad_options (OPT, VAL)
  992.      When called with two arguments, this function allows you set
  993.      options parameters for the function `quad'.  Given one argument,
  994.      `quad_options' returns the value of the corresponding option.  If
  995.      no arguments are supplied, the names of all the available options
  996.      and their current values are displayed.
  997.    Here is an example of using `quad' to integrate the function
  998.        F(X) = X * sin (1/X) * sqrt (abs (1 - X))
  999. from X = 0 to X = 3.
  1000.    This is a fairly difficult integration (plot the function over the
  1001. range of integration to see why).
  1002.    The first step is to define the function:
  1003.      function y = f (x)
  1004.        y = x .* sin (1 ./ x) .* sqrt (abs (1 - x));
  1005.      endfunction
  1006.    Note the use of the `dot' forms of the operators.  This is not
  1007. necessary for the call to `quad', but it makes it much easier to
  1008. generate a set of points for plotting (because it makes it possible to
  1009. call the function with a vector argument to produce a vector result).
  1010.    Then we simply call quad:
  1011.      [v, ier, nfun, err] = quad ("f", 0, 3)
  1012.           => 1.9819
  1013.           => 1
  1014.           => 5061
  1015.           => 1.1522e-07
  1016.    Although `quad' returns a nonzero value for IER, the result is
  1017. reasonably accurate (to see why, examine what happens to the result if
  1018. you move the lower bound to 0.1, then 0.01, then 0.001, etc.).
  1019. File: octave,  Node: Orthogonal Collocation,  Prev: Functions of One Variable,  Up: Quadrature
  1020. Orthogonal Collocation
  1021. ======================
  1022.  - Loadable Function: [R, A, B, Q] = colloc (N, "left", "right")
  1023.      Compute derivative and integral weight matrices for orthogonal
  1024.      collocation using the subroutines given in J. Villadsen and M. L.
  1025.      Michelsen, `Solution of Differential Equation Models by Polynomial
  1026.      Approximation'.
  1027.    Here is an example of using `colloc' to generate weight matrices for
  1028. solving the second order differential equation U' - ALPHA * U" = 0 with
  1029. the boundary conditions U(0) = 0 and U(1) = 1.
  1030.    First, we can generate the weight matrices for N points (including
  1031. the endpoints of the interval), and incorporate the boundary conditions
  1032. in the right hand side (for a specific value of ALPHA).
  1033.      n = 7;
  1034.      alpha = 0.1;
  1035.      [r, a, b] = colloc (n-2, "left", "right");
  1036.      at = a(2:n-1,2:n-1);
  1037.      bt = b(2:n-1,2:n-1);
  1038.      rhs = alpha * b(2:n-1,n) - a(2:n-1,n);
  1039.    Then the solution at the roots R is
  1040.      u = [ 0; (at - alpha * bt) \ rhs; 1]
  1041.           => [ 0.00; 0.004; 0.01 0.00; 0.12; 0.62; 1.00 ]
  1042. File: octave,  Node: Differential Equations,  Next: Optimization,  Prev: Quadrature,  Up: Top
  1043. Differential Equations
  1044. **********************
  1045.    Octave has two built-in functions for solving differential equations.
  1046. Both are based on reliable ODE solvers written in Fortran.
  1047. * Menu:
  1048. * Ordinary Differential Equations::
  1049. * Differential-Algebraic Equations::
  1050. File: octave,  Node: Ordinary Differential Equations,  Next: Differential-Algebraic Equations,  Prev: Differential Equations,  Up: Differential Equations
  1051. Ordinary Differential Equations
  1052. ===============================
  1053.    The function `lsode' can be used Solve ODEs of the form
  1054.      dx
  1055.      -- = f (x, t)
  1056.      dt
  1057. using Hindmarsh's ODE solver LSODE.
  1058.  - Loadable Function:  lsode (FCN, X0, T, T_CRIT)
  1059.      Return a matrix of X as a function of T, given the initial state
  1060.      of the system X0.  Each row in the result matrix corresponds to
  1061.      one of the elements in the vector T.  The first element of T
  1062.      corresponds to the initial state X0, so that the first row of the
  1063.      output is X0.
  1064.      The first argument, FCN, is a string that names the function to
  1065.      call to compute the vector of right hand sides for the set of
  1066.      equations.  It must have the form
  1067.           XDOT = f (X, T)
  1068.      where XDOT and X are vectors and T is a scalar.
  1069.      The fourth argument is optional, and may be used to specify a set
  1070.      of times that the ODE solver should not integrate past.  It is
  1071.      useful for avoiding difficulties with singularities and points
  1072.      where there is a discontinuity in the derivative.
  1073.    Here is an example of solving a set of three differential equations
  1074. using `lsode'.  Given the function
  1075.      function xdot = f (x, t)
  1076.      
  1077.        xdot = zeros (3,1);
  1078.      
  1079.        xdot(1) = 77.27 * (x(2) - x(1)*x(2) + x(1) \
  1080.                  - 8.375e-06*x(1)^2);
  1081.        xdot(2) = (x(3) - x(1)*x(2) - x(2)) / 77.27;
  1082.        xdot(3) = 0.161*(x(1) - x(3));
  1083.      
  1084.      endfunction
  1085. and the initial condition `x0 = [ 4; 1.1; 4 ]', the set of equations
  1086. can be integrated using the command
  1087.      t = linspace (0, 500, 1000);
  1088.      
  1089.      y = lsode ("f", x0, t);
  1090.    If you try this, you will see that the value of the result changes
  1091. dramatically between T = 0 and 5, and again around T = 305.  A more
  1092. efficient set of output points might be
  1093.      t = [0, logspace (-1, log10(303), 150), \
  1094.              logspace (log10(304), log10(500), 150)];
  1095.  - Loadable Function:  lsode_options (OPT, VAL)
  1096.      When called with two arguments, this function allows you set
  1097.      options parameters for the function `lsode'.  Given one argument,
  1098.      `lsode_options' returns the value of the corresponding option.  If
  1099.      no arguments are supplied, the names of all the available options
  1100.      and their current values are displayed.
  1101.    See Alan C. Hindmarsh, `ODEPACK, A Systematized Collection of ODE
  1102. Solvers', in Scientific Computing, R. S. Stepleman, editor, (1983) for
  1103. more information about the inner workings of `lsode'.
  1104.